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ABSTRACT 



This thesis presents a Fortran program that numerically solves the 
steady-state matrix Riccati equation of the quadratic cost optimal 
control problem. Each step of the program is presented, analytically 
and computationally. The check points incorporated in the program and 
the input parameters that can be used to assure a correct solution are 
identified and discussed. Difficulties encountered when verifying the 
program, and the suggested solutions, are also presented. 
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I. INTRODUCTION 



The Fortran computer program presented in this thesis provides a 
relatively quick and reliable means for determining the unique, 
symmetric, steady-state solution P of the nonlinear matrix Riccati 



equation: 


P = 0 = - PA - A'P - C'C + PBR - 1 B 1 P (1) 



occurring in optimal control theory. The remaining equation variables 
are defined in the following statement of the quadratic cost optimal 
control problem. 

Given the linear, time-invariant, completely controllable system 
defined by the state equations: 



where: 


x (t) = Ax(t) + Bu(t) (2) 

x (0) = x Q (3) 

x(t) is an n x 1 state vector 

u(t) is an m x 1 unconstrained control vector 

A is an n x n matrix 

B is an n x m matrix. 



determine the control vector u*(*) which minimizes the quadratic cost 
functional : 

J(x ;u(- ) ) = J [ x ' (t)C 'Cx(t) + u ' (t)Ru(t) ] dt (4) 

o ' 7 



where: 


C is an m x n matrix 

R is an m x m symmetric, positive definite matrix 



and the matrix pair (A,C) is completely observable. 
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It has been shown [Ref. 1] that u*(-) is given by the linear feed- 
back law: 

u*(t) = -R~^B'Px(t) = -L*x(t) (5) 

where P is the unique solution to matrix equation (1). 

For the system of equation (2) to be completely controllable, the 
n x mn augmented matrix G, 

G = (B | AB | A 2 B | ••• | A n_1 B ), (6) 

must contain n linearly independent column vectors, or equivalently, 
the rank of G must be n. For the matrix pair (A,C) to be completely 
observable, the n x mn augmented matrix H, 



H = ( C' | A 1 C * | (A') 2 C | ••• | (A’) n " 1 C' ), 

must contain n linearly independent column vectors or the rank of H 
must be n. 

When u*( - ) is given by equation (5), equation (2) can be rewritten 
as: 

x(t) = Ax(t) - BL*x(t) = (A - BL*)x(t) (7) 



which, using equation (3), has the general solution: 

x(t) - e< A - BL *)‘ . . 



( 8 ) 



Because the system is completely controllable, as time approaches 
infinity x(t) must remain bounded. To satisfy this requirement, the 
matrix (A - BL*) must be a stable matrix or, alternatively, all the 
eigenvalves of the matrix (A - BL*) must have negative real parts. As 
a consequence of this: 

1 im 



t x(t) = x„ e^ A - BL *> = 0 



(9) 
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To evaluate the quadratic cost associated with the optimal con- 



trol u*(t), substitute equation (5) into equation (4) and write: 

/ oo 

[x ' (t)C'Cx(t) + x'(t)L*'RL*x(t)] dt. (10) 

Substituting equation (8) for x(t), expanding L*'RL*, and factoring out 
the common terms: 



J* = x' 
o 



r" e (A - BL*)'t [c , c + p BR -l B , p] e (A - BL*)t dt ‘ 

0 J 



x 0 . (ID 



From equation (1) C'C = -PA -A'P + PBR B'P or: 



”1 n I | 



C'C + PBR -1 B 1 P = -P (A - BL*) - (A - BL*)'P . 



(12) 



Substituting equation (12) into equation (11) and integrating by parts 
the first term of the resulting integral expression, we get: 



f e^ A ' BL *)' t [— P (A - BL*)] e^ A ' BL *^ dt 

o J 



e (A - BL*) ' t [p] e (A - BL*)t 



t=0 



/'VBL*)'e< fl - BL *>' t [P]e‘ A - BL *) t dt. 



(13) 



Using equation (9) to evaluate the upper limit of the first term, the 
first term reduces to P. Recognizing: 

r 

Ze Zt = l l JT ZV = e Zt Z (14) 

i=0 1 * 

we see that the second term in equation (13) will cancel the second 
term in equation (11) when equation (12) is substituted. Thus: 



J* = 



x'Px 

o 



o * 



(15) 
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II. ANALYTICAL METHODS OF SOLUTION 



A. ANALYTICAL METHOD OF KLEINMAN 

The main computer program of this thesis is based upon a method for 
solving the steady-state Riccati equation published by David L. Kleinman 
[Ref. 2] in 1968. The iteration scheme for solving euqation (1) is 
presented in the following theorem by Kleinman. 

1 . Kleinman 's Theorem 

Let V k , k = 0, 1 , 2 , . . . be the n x n (unique) positive definite 
matrix solution of the linear algebraic matrix equation: 

0 = A k V k + V k A k + C ' C + L k RL k (16) 

where, recursively. 



* R_1 B ' v k-1 


k = 1,2,3,... 


. A - BL k 


k = 0,1,2,... 



and where L„ is chosen such that the matrix A„ = A - BL„ has eigenvalues 
o o o 3 

with negative real parts. 

Then: 1) P < V k+] < V k < • • • k = 0,1,2,... 

pi 1 im V = P 
L) k-*» v k K * 

The notation X > Y [X > Y] means that the matrix X - Y is 
positive [semi] definite. 

2. The Cost Matrix 

The proof of this theorem requires the introduction and defini- 
tion of a cost matrix V^. Assume that u k (x(t)) = - Lx(t) is an arbitrary 
feedback law, with feedback gains of L, and u k (x(t)) is applied to the 
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system of equation (2). Following a development similar to equations 
(7), (8) and (15), the resulting quadratic cost function is: 



J(x 0 ;u L (x(t))) = x^V L x 0 

where is the cost matrix associated with the feedback gains L and is 
given by: 

V, = f e^ A " BL ^ ,t: [C'C + L'RL] e^ A " BL ^ dt . (17) 

o J 

is bounded if and only if the closed-loop system control matrix (A - BL) 
is stable. If is bounded, it becomes the unique (positive definite) 
solution of the linear matrix equation: 

0 = (A - BL)'V l + V l (A - BL) + C'C + L'RL . (18) 

Examining the first term of equation (18) and substituting 
equation (17) for V^, we can verify this relationship: 

(A - BL) ' V = / (A - BL)' e^ A ' BL ^ ^CC'C + L'RL] e^ A " BL ^ dt. 

L o J 

Integrating by parts and using equation (14) we have: 

co 

(A - BL ) ' V 1 = e^ A ‘ [C'C + L'RL] e^ A ‘ BL ^ t=Q 

f e (A - BL)'tj- c , c + L , RL -j e (A - BL)t ^ _ BL j dt> 

o J 

Using equation (9) to evaluate the upper limit of the first term, the 
first term reduces to C'C + L'RL and we have, 

(A - BL)'V l =-C'C - L'RL - V l (A - BL) , 

the desired result. 
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Recalling the result of equation (15) we see that for the 
optimal control of equation (5) we have: 

V L * = P . (19) 



If L-j and L 2 are the gain matrices associated with the cost 
matrices V-| and V 2 , it can be shown [Ref. 3] that: 



V 



1 




(A - BL 2 ) 1 1 



[(L 1 - L 2 )'R(L 1 




or: 



- (^ - l 2 ) , (b , v 1 - rl 2 ) 

- (B'V 1 - RL 2 ) , (L 1 - L 2 )] e^ A ‘ BL 2^ dt 




V 2 = f e (A ' BL l },t [(L 1 - L 2 )'R(L 1 - 

- (L 1 - L 2 )'(B'V 2 - RL 2 ) 

- (B'V 2 - RL 2 ) ' (L 1 - L 2 )3 e (A ' BL l }t 




dt. 



( 20 ) 



( 21 ) 



If either matrix (A - BL 2 ) or matrix (A - BL-j) is unstable, V-j 
or V 2 , respectively, will be unbounded and care must be exercised in 
using the above formulas. 

3. Proof of Kleinman's Theorem 

1) P < V R+1 < V R s ••• k = 0,1,2, •••. Let V Q be the cost 
matrix for L o , the initial guess that yields a stable closed-loop 
system control matrix, and let be the cost matrix for L-| = R~^B'V 0 . 
Substituting V Q and V-| into equation (20) and noting: 



e A 'k t = (e A k t ) ' 



R = Y'Y, where Y is unique 



B'V 0 - RL 1 = 0 
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we have: 



v o - V 



! = / e A l t [(L o - L-, ) 'R(L o - L-j )] e A l t dt 



A n t 



Define Z(t) = Y(L Q - L-|)e 1 . It has been shown [Ref. 4] that for 
Z(t) a real matrix: 

Z'(t)Z(t) * 0 for all t i 0. 

/ co 

Z* (t)Z(t) dt > 0 or: 

V v i • 



( 22 ) 



Let be the cost matrix associated with L*, use equation (19) 
and substitute into equation (21). Noting that: 



B'P - RL* = 0 



we have: 



- P = / e A l' t [(L 1 - L*)'R(L 1 - L*)] e A l t dt. 



Thi 



At t 



is time define Z(t) = Y(L^ - L*)e 1 and we have: 



V, - P = / Z'(t)Z(t) dt > 0 or: 

O'' 



V ] ^ P . (23) 



Combining the results of equations (22) and (23) we get: 



V o * V 1 * P 



meaning that V-| is bounded above by P and below by V Q . Thus, A-| is a 
stable matrix and satisfies equation (16) with k = 1. Similar 
arguments can be made for k = 2,3,4,... yielding: 
v 0 2 V, > V 2 > ••• V k > V k+1 > P . 

the desired result. 
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2) ^ V k = P. The V k = V^exists by a theorem on monotonic 

convergence of positive operators [Ref. 5], Thus, in the limit = 
V k+1 and A k = A - BL k = A - BR _1 B'V k = A - BR _1 B'V k+1 . Since A£ = 

A' - V k BR"^B' equation (16) becomes, in the limit as k approaches 
infinity, 

0 = A'V k - V k BR _1 B'V k + V k A - V k BR _1 B'V k + C'C + V k BR' 1 B , V k 
which is equation (1), the desired result. 

B. BASS'S METHOD OF DETERMINING A STABLE INITIAL GUESS 

Kleinman's method requires an initial guess of the feedback gain 
that yields a stable closed-loop system control matrix. Kleinman 
remarks [Ref. 2] that, if the system of equation (2) is completely 
controllable, then the desired initial guess will always exist. A 
method that could be programmed for a general, controllable system to 
yield this correct initial guess was sought. 

Bass [Ref. 6] presented, but did not publish, a general method for 
determining a stable initial guess for a completely controllable system 
in 1961. The method was published in a paper by Wonham and Cashman 
[Ref. 7] and a proof can be found in a subsequent paper by Bass [Ref. 8]. 

Given the control labe matrix pair (A,B), the matrix: 

A rt = A - BL 
o o 

will be stable when L Q = B'X”\ where X is the (unique) positive definite 
solution to the linear matrix equation: 

- (A + $ 1 ) X - X(A + 31 )' + 2BB' = 0 (24) 

where B is defined by: 

B > I I A | | 
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where ||*|| is the Euclidean norm. According to Wonham [Ref. 7] the 
results are also valid if: 



0 att- M AX r , a 

2 > . I I a.. 

J i=l 1J 



where a., are the elements of A. 

* J 



C. SOLVING THE MATRIX EQUATION Y'X + XY = Z 

Equations (16) and (24) are in the form of the Lyapunov equation: 

Y'X + XY = Z (25) 

The methods found in the literature for solving this equation fall into 
two general categories: series solutions and simultaneous linear 

equation solutions. 

In the series solutions, X is found from the sum of a converging 
matrix series, i.e.. Ref. 9. For the series to converge the matrix Y 
must be a stable matrix. In solving equation (16) this condition is 
met; A^ is stable by the definition of a bounded cost matrix. In 
general, however, the matrix (A + Bl)' of equation (24) is not stable, 
and the series method fails to solve equation (24) properly. 

Since the unique solution, X, is symmetric, equation (25) represents 
n(n + l)/2 unknowns. The second category of solutions expands equation 
(25) into a set of n(n + l)/2 simultaneous linear equations. An 
economical way of recursively expanding equation (25) was given by 
Bingulac [Ref. 10], using an integer coefficient matrix to expand the 
equation. This is the method used to expand equations (16) and (24) 
to the form Ax = B, which can be solved by a variety of simultaneous 
equation solvers. 
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D. ALTERNATE ANALYTICAL METHODS 



There are several alternate methods found in the literature for 
solving the steady-state matrix Riccati equation (1). 

A method developed by Bass [Ref. 11] obtains the solution from a 
2n-dimensional Hamiltonian, H, and the terms of the polynomial expansion 
of the stable roots of H. This scheme is considered too sensitive to 
finite numerical computations to be of practical use. 

MacFarlane [Ref. 12] shows that the solution can be obtained from 
the eigenvectors corresponding to the unstable eigenvalues of a 
similar Hamiltonian. The scheme requires that the system have distinct 
eigenvalues and that the corresponding eigenvectors be determined (this 
is difficult for high order systems). 

Blackburn [Ref. 13] programmed a method based on a Newton-Raphson 
iterative solution for simultaneous nonlinear equations. This scheme 
requires an initial guess that yields a stable closed-loop system con- 
trol matrix. The user must determine this initial guess so that it is 
close enough to the final solution for the local Newton-Raphson method 
to converge. 

A fourth method integrates the full, nonsteady-state Riccati 
equation (1) backwards in time, from a set of zero initial conditions, 
until each element of the P matrix reaches a satisfactory, small value. 
For systems of even moderate order, this method is prohibitive with 
respect to computation time. 
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III. COMPUTATIONAL METHOD OF SOLUTION 



A. SUBROUTINE RICATS 

1 . General 

Subroutine RICATS is the subroutine that the user will call to 
solve equation (1). The program language is FORTRAN IV for the Operat- 
ing System / 360 which is compatible with, and encompasses USASA 
FORTRAN. The calling arguments, in the required sequence, are explained 
in the comment cards at the beginning of the subroutine. It should be 
noted that IA = n and JB = m. Basically the subroutine iterates on 
equation (16), using equation (24) to calculate the initial guess, until 
the solution converges. 

For first order systems (n = 1), the nonlinear equation (1) 
becomes a quadratic equation. For these systems, subroutine RICATS 
solves the resulting quadratic and returns the largest root in the 
first element of the P matrix . 

2. The System Controllability Check 

An analytical requirement of Bass's method of determining the 
initial guess is that the matrix pair (A,B) be completely controllable. 
To check this requirement the augmented matrix G of equation (6) is 
formed and subroutine GMRANK is called to determine the rank of G. If 
the rank equals n, execution continues, otherwise subroutine RICATS 
returns with IER = 2. 

3. The Initial Guess Stability Check 

Kleinman's iteration scheme requires that the eigenvalues of the 
closed-loop system control matrix of the initial guess have negative 
real parts in order for the initial cost matrix V Q to be bounded. 
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In subroutine RICATS, when the stability check is requested, the matrix 
A - BL q is copied into a work matrix and then passed to the IBM/SSP 
subroutines HSBG and ATEIG. The eigenvalues computed by ATEIG are then 
individually tested to be algebraically less than minus the absolute 
value of EIGMAX, where EIGMAX is set by the user. If any eigenvalue 
fails the check RICATS returns with IER = 3. 

4. The Solution Positive Definite Check 

Kleinman's iteration scheme is supposed to converge to a positive 
definite solution. The iterations can converge to a nonpositive definite 
solution if all of Kleinman's theorem requirements are not numerically 
met. Therefore, for the user's convenience, a check to verify that the 
solution is positive definite is included. In subroutine RICATS, when 
the positive definite check is requested, the solution is factored by 
the Cholesky square-root method. A nonpositive term on the diagonal 
of the factorized matrix constitutes a nonpositive definite solution 
and RICATS returns with IER = 4 + KL, where KL = 0 is for a converged 
solution and KL = 1 is for NTRY iterations without convergence. 

5. The Steps of Subroutine RICATS 

Calling subroutine RICATS causes the following sequential 
steps to be executed. 

1. Check input parameters IA, JB, IER, NTRY to insure they are 
within proper bounds. If check fails IER = 6, NN = 0 and RICATS returns. 

2. If the user requests, check the controllability of the 
inputed system. If check fails IER = 2, NN = 0 and RICATS returns. 

3. Set NN = 0. 

4. If this is a first order system (n = 1), go to step 30. 

5. Form E = - BR _ 1 B ' . 
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6. Form F = 2BB ' . 

7. Form P = (A + gl)' using equation (24) to define e where 
the "greater than" magnitude is provided by the variable FIX or: 



8. Solve (A + eI)X + X(A + el)' = 2BB * . 

9. If subroutine SIMQ, through subroutine MLIAPS, returns a 
nonzero error flag, NN = 1 . 



1 1 . Form P = A - BL Q . 

12. If the user requests, check the stability of the system 
matrix of the initial guess. If the check fails, IER = 3 and RICATS 
returns. 

13. Form F = - Q - L'RL . 

^ oo 



15. If subroutine SIMQ, through subroutine MLIAPS, returns a 
nonzero error flag, NN = NN + 1 . 

16. Set k = 1, KL = 0. 

17. Form P = EV k . 

18. Form F = - Q + V^P. 

19. Form P = A + P. 

20. Solve (A + EV k )'V k+1 + V k+] (A + EV k ) = - Q + V k EV R . 

21. If subroutine SIMQ, through subroutine MLIAPS, returns a 
nonzero error flag, NN = NN + 1 . If NN > (n + l)/2, go to step 29. 

22. Check each element of V k+ ^ by ABS ( (v k - v k+ -|)/v k ) 

< TOLER. If all elements of V k+ ^ pass this test, go to step 27. 

23. Form » k = V k+ , . 




n 



10. FormL 0 =BX _1 . 



14. Solve (A - BL 0 )'V 1 + V-, (A - BL Q ) = - Q - L'RL^ 
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is a positive 



24. If k > NTRY; go to step 26. 

25. Set k = k+1 ; go to step 17. 

26. Set KL = 1. 

27. If the user requests, check to see if 
definite matrix. If check fails IER = 4 + KL and RICATS returns. 

28. Set IER = KL and RICATS returns. 

29. Set IER = 7 and RICATS returns. 

30. If the user requests the controllability check and B(l) is 
zero, IER =2, NN = 0 and RICATS returns. 

31 . Set P = AR/BB + ^(AR/BB) 2 + C'CR/BB 

32. If the user requests the positive definite check and P s 
0.1E-35, IER = 4, NN = 0 and RICATS returns. 

33. Set IER = 0 and RICATS returns. 

B. SUBROUTINE MLIAPS 

Subroutine MLIAPS is an auxiliary routine used by subroutine RICATS. 
The subroutine expands the equation in steps 8, 14 and 20 of subroutine 
RICATS into n(n + l)/2 simultaneous linear equations of the form Ax = B. 
The method used was presented by Bingulac [Ref. 10] in 1970. These 
n(n = l)/2 simultaneous linear equations are then solved by subroutine 
SIMQ. Upon return from subroutine SIMQ, MLIAPS immediately returns to 
subroutine RICATS and any error codes returned by SIMQ are passed, 
unchanged, to RICATS. 

C. SUBROUTINE GMRANK 

Subroutine GMRANK is a general subroutine for determining the mini- 
mum row or column rank of an arbitary real matrix. The method used is 
simple row and column interchanges for maximum pivoting, with successive 
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reduction on the remaining matrix elements [Ref. 14]. When the absolute 
value of a pivot element is less than 0.1E-35, or when the final pivot 
element has been determined, the subroutine returns the integer K, 
where K is the number of successfully determined nonzero pivot elements 
or, equivalently, the rank of the inputed matrix. 
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IV. OPERATIONAL ASPECTS OF USING SUBROUTINE RICATS 



A. CONVERGENCE PROPERTIES 

Kleinman suggests that equation (16) will converge to a satisfactory 
solution within seven to ten iterations, using an initial guess generated 
by hand or from a previous solution. Subroutine RICATS, using a correct 
set of input parameters and Bass's method of generating an initial guess, 
converged within forty iterations for every system tested. 

To test subroutine RICATS, over seventy-five systems were run and a 
satisfactory solution was returned for each one. As a test of how 
accurate the solutions were, the average absolute value of the matrix P 
of equation (1) was determined. For the satisfactory solutions the 
average values were of the order of 10”^. The standard input parameters 
were: 

NTRY = 50 

TOLER = 0.1E-03 
FIX = 1.1 

EIGMAX - 0.001 
IER = 0 

and the usual error flag returned was IER = 0 (see the discussion of 
the parameter IER for high order systems). 

B. INPUT PARAMETERS 

Through its input parameters subroutine RICATS was written to be as 
flexible as possible for the user. For any system, the five parameters 
NTRY, TOLER, FIX, EIGMAX, IER influence the final returned solution. 

It is hoped that the user can tailor these input parameters to meet 
his particular needs. 

In the following discussion, "high" or "higher order systems," 
refer to systems of the eighth order and above (n > 8). 
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1 . NTRY 



NTRY is the maximum number of iterations the user desires to 
be attempted, before returning to the user's program without a converged 
solution. A recommended value is fifty and the user is usually assured 
that a solution that is going to converge, will converge within fifty 
iterations. It should be noted that some initial guesses generated may 
yield marginally stable system control matrices, and these solutions will 
require a larger number of iterations to converge. From experience, one 
hundred and fifty iterations was considered to be a sufficient practical 
upper 1 imit. 

2. TOLER 

TOLER is the maximum percentage difference, between each element 
of the solution, on successive iterations for the convergence criterion 
to be satisfied. The accuracy of single-precision computations is about 
five significant digits, so decreasing TOLER beyond 0.0001 will not 
result in an appreciably more accurate solution. Decreasing the magni- 
tude of TOLER (up to this limit) will tend to increase the accuracy 

and the number of iterations required for the desired solution. 

/ 

3. FIX 

FIX is the constant used to insure the magnitude of g in equa- 
tion (24) is greater than the Euclidian norm of the A matrix. 
Analytically then, FIX should be greater than one and the suggested 
value of 1.1 worked well for the majority of systems. However, for 
some systems equation (24) generates a singular initial guess, evidenced 
by an underflow error message from subroutine MINV. These singular 
initial guesses can still yield quite satisfactory solutions within 
five or six iterations. However, the user can, by varying FIX through 
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the suggested range, remove this underflow from his execution. From 
practical experience, the range of FIX is from 0.1 up to about three 
or four. 

4. EI6MAX 

Once the initial guess L Q is calculated, the user has the option 
through IER of verifying that the associated control matrix is stable. 

If the user asks for the stability check, the eigenvalues of the control 
matrix are determined. Each must have its real part more negative than 
minus the absolute value of EIGMAX. From the discussion on cost matrices, 
theoretically EIGMAX could be set to zero. Numerically though, an 
eigenvalue that is very close to zero can induce numerical instability 
in the iteration scheme. From experience, the suggested value is 0.001. 
The user can also use FIX to modify the control matrix of the initial 
guess in an attempt to make the real parts of the eigenvalues negative 
enough to pass the stability check. 

5. IER 

The most informative parameter of the group is IER. On return 
from RICATS, IER can inform the user of the validity of the solution. 

When speicfying IER as in input parameter, the user can, by bypassing 
various combinations of the three checks available in RICATS, overcome 
some system deficiencies, save execution time, or obviate decks for 
the subroutines GMRANK, ATEIG, HSBG (see note 6 in comment cards at the 
beginning of subroutine RICATS). By bypassing any or all of the three 
checks (if the user is fairly certain that the circumvented checks 
would have been passed ) the user can save execution time, although 
the savings are neglibile except for high order systems. 
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Care must be exercised when not checking controllability and 
or stability. For systems that are uncontrollable, neglecting the 
controllability check may lead to a convergent solution, if the stability 
check were passed. Such a solution could be positive definite and yield 
a stable final-solution control matrix, but the user should keep in 
mind that there may be modes of the system that cannot be controlled. 

For the case of unstable control matrices corresponding to the initial 
guess, the stability check should be bypassed only after FIX has been 
varied through the suggested range of values without success. The 
danger in attempting to generate a solution for a system that would fail 
the stability check is that successive iterations are no longer bounded 
by a lower positive [semi] definite iteration. The iterations will 
probable not converge to a satisfactory answer. This is the most 
dangerous of the three checks to bypass, by far. 

The positive definite check of the solution is included as a 
convenience for the user, to verify that the final solution is indeed 
positive definite. A note of caution should be sounded when solving 
high order systems. The positive definite check is subject to numerical 
problems when evaluating the high order solutions, and thus a solution 
will be flagged as not being positive definite, when for all practical 
purposes it is. Refinements in the solution, from introducing iterative 
routines for solving equation (16) (see subroutine SIMQ below), have 
shown that the difference between passing and failing the check can 
depend solely on numerically insignificant digits in a few elements of 
the solution. Therefore, to give confidence to a solution that has been 
flagged as not being positive definite, the user can: (1) determine 

the eigenvalues of the final solution, closed-loop system control 
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matrix and check for all negative real parts; or (2) calculate the value 
of the optimal quadratic cost function for an arbitrary set of initial 
conditions using equation (15), and check for a positive, finite result. 
The positive definite check will yield the same result as the IBM/SSP 
subroutine MFSD. 

C. HIGH ORDER SYSTEMS 

As the order of the system increases, the problems due to finite 
numerical calculations increase considerably. One means of trying to 
maintain sufficient accuracy is to have a double-precision version of 
subroutine RICATS. Since this would nearly double the storage require- 
ments (prohibitive for systems of order higher than twenty) it was felt 
that this was not a feasible means of achieving the goal. The suggested 
procedure for systems of higher order is for the user to introduce his 
own dummy subroutines MINV and SIMQ. Care must be exercised when writing 
your own subroutines. The variables passed to and returned from your 
subroutine must correspond exactly to the variables as handled by the 
subroutine you are replacing. 

1 . Subroutine MINV 

Using a common statement from the user's main program to provide 
the additional storage required, the user can write a routine to convert 
the matrix to be inverted to double-precision storage, invert the matrix 
in double-precision, then convert the inverted matrix back to the 
single-precision mode. When this is done, the inverted matrix is 
passed back to the subroutine RICATS as if the IBM/SSP subroutine MINV 
had done the inverting. A logical way to invert the matrix in the double- 
precision mode is to use the IBM/SSP subroutine DMINV. It might appear 
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that this type of routine would make no noticeable change in the 
execution of subroutine RICATS. However, for some systems tested, 
particularly those with singular initial guesses (FIX = 1.1), the number 
of iterations required for solution was halved. 

2. Subroutine SIMQ 

Again using a common statement, recognizing that work areas can 
be declared single-precision in one subroutine and double-precision in 
another, the user can write his own simultaneous linear-equation-solving 
routine. The IBM/SSP subroutine SIMQ has been found to yield satis- 
factory results for a maximum of about forty-five equations (n = 9). 

For systems of higher order, the solutions did converge, but they were 
usually accompanied by an error flag indicating that they were not 
positive definite solutions. However, iterative refinement of the 
solution to equation (16) at each iteration can yield error-free 
results. The user can either write his own routine for the iterative 
refinement, or pass the set of equations to the IBM/SSP subroutines 
FACTR and RSLMC. It is suggested that the common statement also contain 
a relative tolerance parameter, not necessarily the same as TOLER 
inputed to RICATS, in either name or magnitude. As was previously 
mentioned, this subroutine technique was used to remove error flags 
from the solutions returned by subroutine RICATS. 
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oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooonoooo 



SUBROUTINE RICATS 



PURPOSE 

TO ITERATIVELY DETERMINE THE UNIQUE, SYMMETRIC, 
POSITIVE-DEFINITE SOLUTION P TO THE NONLINEAR, 
STEADY-STATE MATRIX RICCATI EQUATION 
. -I 

P = 0 = - P* A - A'*P - C'*C + P*B*R *B'*P 
OF OPTIMAL CONTROL THEORY. SEE NOTE 10. 



USAGE 

CALL 

1 



RICATS (A,6,Q,R,P,D,E,F,V,L,LI,MI,IA,JB,KQ, 
NTRY, TOLER, F IX, E I GMAX , I E R , NN ) 



DESCRIPTI ON 
A 

B 

Q 



D 

E 

F 

V 

L 

LI 
MI 
I A 



JB 

KQ 



IER 



IA INPUT MATRIX. SEE NOTE 



OF PARAMETERS 
GENERAL IA BY 

3. 

GENERAL IA BY JB INPUT MATRIX. 

LOWER TRIANGULAR (OR UPPER TRIANGULAR, 
SEE KQ) PART OF A SYMMETRIC, PCSITVE 
SEMIDEFINITE IA BY IA INPUT MATRIX. 

Q = C ' C . 

GENERAL PART OF A SYMMETRIC, POSITIVE 
DEFINITE JB BY JB INPUT MATRIX. 

IA BY IA WORK MATRIX. GN RETURN P CON- 
THE GENERAL FORM OF THE SOLUTION. 



TAINS 
MM BY 



MM 
MZ 
MM 
I A 
I A 
I A 
ROW 



BY 
BY 
BY 
BY 
BY 
BY 
AND 



MM WORK MATRIX. 
1 WORK VECTOR. 
1 WORK VECTOR. 
1 WORK VECTOR. 
IA INTEGER WORK 
1 INTEGER WORK 
1 INTEGER WORK 



SEE NOTE 1, 

SEE NOTE 1 

MATRIX. 
VECTOR. 
VECTOR. 



COLUMN DIMENSION OF MATRICES A, 



NTRY 

TOLER - 



FIX 



El GMAX - 



Q AND RETURNED SOLUTION P. ROW 
DIMENSION OF MATRIX B. 

COLUMN DIMENSION OF MATRIX S. ROW AND 
COLUMN DIMENSION OF MATRIX R. JB < = I A. 
INTEGER INPUT CONSTANT. 

= 0 - Q MATRIX IS STORED COLUMNWISE IN 
LOWER TRIANGULAR FORM. 

= 1 - Q MATRIX IS STORED COLUMNWISE IN 
UPPER TRIANGULAR FORM (SEE 
I3M/SSP SUBROUTINF MSTR). 

MAXIMUM NO. OF ITERATIONS TO BE TRIED. 

2 > NTRY > 151. 

INPUT CONSTANT WHICH IS USED AS A 
RELATIVE TOLERANCE FOR TEST OF 
CONVERGENCE. 

INPUT CONSTANT. THEORETICALLY SHOULD BE 
>= 1.0. SUGGESTED VALUE IS 1.1. SEE 
NOTE 7. 

MINIMUM ACCEPTABLE MAGNITUDE OF THE 
NEGATIVE REAL PARTS OF THE EIGENVALUES 
OF THE CONTROL MATRIX ( A - B*V(0) ) OF 

THE INITIAL GUESS V(0). SUGGESTED VALUE 
OF EIGMAX IS OF ORDER 0.001 . SEE NOTE 

5. EIGMAX USED ONLY FOR IER = 0,1, 4, 5. 
INTEGER INPUT/OUTPUT PARAMETER. 

THE FOLLOWING NOTATION IS USED BELOW : 
FOR NUMERICAL CONTROLLABILITY 
OF THE MATRIX PAIR (A, 8). 

SEE NOTES 4,6. 

FOR A STABLE CONTROL MATRIX FOR 
THE INITIAL GUESS. SEE NOTES 5,6. 
FOR A POSITIVE-DEFINITE SOLUTION. 
SEE NOTE S. 



C . A , B . 



S.C.M. - 
P.D.S. - 
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NN 



INPUT - THE VALUE OF IER DETERMINES WHICH 
OF THE 3 ABOVE PROPERTIES RICATS 
IS TO CHECK DURING EXECUTION. SEE 



— 


0 - 


NOTE 9. 
CHECK C . A , B . 


; S.C.M. 


; P.D.S. 




1 - 


CHECK C .A, B. 


; S.C.M. 






2 - 


CHECK C . A , B . 




P.D.S. 


— 


3 - 

4 - 


CHECK C . A , B . 
CHECK 


S.C.M. 


5 P . D . S . 


= 


5 - 


CHECK 


S.C.M. 




= 


6 - 


CHECK 




P.D.S. 


= 


7 - 


MAKE NONE OF 


THE ABOVE 


3 CHECKS 



OUTPUT - ON RETURN IER SERVES AS AN ERROR 
FLAG. 

= 0 - ALL CHECKS REQUESTED WERE PASSED, 
SOLUTION CONVERGED WITHIN NTRY 
ITERATIONS. 

= 1 - ALL CHECKS REQUESTED WERE PASSED, 
SOLUTIONS DID NOT CONVERGE WITHIN 
NTRY ITERATIONS. 

= 2 - C . A , B . CHECK FAILED. NO SOLUTION. 

= 3 - S.C.M. CHECK FAILED. NO SOLUTION. 

= 4 - P.D.S. CHECK FAILED. SOLUTION DID 

CONVERGE WITHIN NTRY ITERATIONS. 

= 5 - P.D.S. CHECK FAILED. SOLUTION DID 
NOT CONVERGE BY NTRY ITERATIONS. 

= 6 - IMPROPER PARAMETERS PASSED TO 

RICATS. SEE NOTE 2. NO SOLUTION. 

= 7 - MORE THAN (IA+D/2 ERROR FLAGS WERE 
RETURNED BY SUBROUTINE SIMQ INDI- 
CATING SINGULAR INPUTS TO THAT SUB- 
ROUTINE. INCORRECT SOLUTION. SEE 
NOT F 7. 

- INTEGER OUTPUT PARAMETER. 

FOR RICATS RETURN WITH IER .NE. 2, NN 
IS THE NO. OF ERROR FLAGS RETURNED 
BY SUBROUTINE SIMQ. 

FOR RICATS RETURN WITH IER .EQ. 2, NN 
IS THE RANK OF THE AUGMENTED 
CONTROLLABILITY MATRIX. SEE NOTE 4. 



NOTES 

1. MM = I A* (IA+D/2 : MZ = MAXO( MM,IA*JB ). 

2. THE INPUT PARAMETERS I A, IB ,1 ER,NTRY MUST 
MEET THE FOLLOWING REQUIREMENTS: 

0 < JB <= IA ; -1 < IER < 8 ; 2 < NTRY < 151 

3. INPUT MATRICES A , B , Q , R AND OUTPUT MATRIX P 
ARE ASSUMED STORED IN THEIR RESPECTIVE MODES 
(GENERAL OR TRIANGULAR) COLUMNWISE IN IBM/SSP 
COMPRESSED FORM. MATRICES A , B , Q , R ARE RETURNED 
UNCHANGED. 

4. THE METHOD OF OBTAINING THE INITIAL GUESS 
REQUIRES THAT THE MATRIX PAIR ( A , B ) IS CONTROLL- 
ABLE (I.E. THE AUGMENTED MATRIX 

2 I A-l 

( B| A-BlA * B | ... | A *B ) 

HAS RANK IA). THE USER CAN HAVE THIS PROPERTY 
CHECKED ( IER=0, 1 , 2 , 3 ) AND IF THE CHECK FAILS, 
RICATS IMMEDIATELY RETURNS WITH IER = 2. 

5. THE ITERATION SCHEME REQUIRES THE CONTROL 
MATRIX OF THE INITIAL GUESS BE STABLE. THE USER 
CAN HAVE THIS PROPERTY CHECKED ( I ER=0, 1,4,5) . 

THE REAL PARTS OF THE EIGENVALUES OF THIS MATRIX 
ARE TESTED TO BE <= -ABS ( E IG MAX ) . IF ANY CHECK 
FAILS RICATS IMMEDIATELY RETURNS WITH IER = 3. 
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6. IF THE USER DOES NOT HAVE C.A,B. CHECKED, 
SUBROUTINE GMRANK IS NOT NEEDED. IF THE USER 
DOES NOT HAVE S.C.M. CHECKED, SUBROUTINES ATEIG 
8 HS8G ARE NOT NEEDED. THE USER CAN OBVIATE UN- 
NECESSARY DECKS BY SUPPLYING HIS OWN DUMMY SUB- 
ROUTINES. AS AN EXAMPLE, GMRANK IS SHOWN BELOW. 

SUBROUTINE GMRANK ( D , I A , KL , TOL ER , K ) 

IA = IA 

RETURN 

END 

THE REQUIRED HEADINGS FOR ATEIG AND HSBG ARE: 
SUBROUTINE ATEIG ( I A , D , V , F , LL , I A ) 
SUBROUTINE HSBG (IA,D,IA) . 

7. CHANGING FIX WILL CHANGE THE INITIAL GUESS. 
VARYING FIX FROM 0.1 TO 4.0 HAS BEEN FOUND TO 
YIELD SATISFACTORY INITIAL GUESSES, THUS REMOVING 
ERROR FLAGS I ER = 3 AND 7, IN MOST CASES. 

8. THE ITERATION SCHEME SHOULD CONVERGE TO A 
POSITIVE DEFINITE SOLUTION. THE USER CAN HAVE 
THIS PROPERTY CHECKED ( I E R=0 , 2 ,4 , 6 ) • FOR THIS 
CHECK THE SOLUTION IS FACTORED BY THE CHOLESKY 
SQUARE-ROOT METHOD. A NONPOSITIVE DIAGONAL ELE- 
MENT OF THE FACTORED MATRIX CONSTITUTES A NON- 
POSITIVE DEFINITE SOLUTION, RICATS RETURNS WITH 

I ER = 4 OR 5. HIGH ORDER ( I A>7 ) SOLUTIONS WITH 
ERROR' FLAGS I ER = 4 OR 5 MATHEMATICALLY MAY BE 
POSITIVE DEFINITE. TO CHECK SOLUTIONS; USE PRIOR 

-1 

EXPERIENCE; FIND EIGENVALUES OF A - B*R *B'*P, 
ALL SHOULD HAVE NEGATIVE REAL PARTS;OR CALCULATE 
Z'*P*Z FOR ARBITRARY IA BY 1 VECTOR Z, RESULTS 
SHOULD BE POSITIVE AND FINITE. 

9. IF ANY OF THE 3 CHECKS FAIL (SEE NOTES 4,5, 
7,8) THE USER SHOULD BYPASS THE FAILED CHECK AND 
ATTEMPT TO FORCE A SOLUTION. BEFORE BYPASSING 
THE S.C.M. CHECK, FIX SHOULD BE VARIED FROM 0.1 
TO 4 AS SUGGESTED IN NOTE 7. 

10. THE STEADY-STATE MATRIX RICCATI EQN . ARISES 
FROM A SYSTEM DEFINED BY THE DIFFERENTIAL EQN. 

X ( T ) = A*X(T) + 6*U(T) , WHERE U(T) IS CHOSEN 
TO MINIMIZE THE QUADRATIC COST FUNCTIONAL J; 

J = INTEGRAL FROM ZERO TIME TO INFINITE TIME OF 
( X* (T)^C' *C*X (T ) + U'(T)*R*Um ) DT . U(CPTIMAL) 

-1 

IS GIVEN BY U(OPTIMAL) = - R *B • *P*X ( T ) . 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED. 

MINV ATEIG SEE NOTE 6. 

SIMQ HSBG SEE NOTE 6. 

GMRANK SEE NOTE 6. 
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SUBROUTINE RICATS { A ,B , Q ,R , P , D , E , F , V , L , L L , M I , I A , J B , KQ, 
1 NTRY, TOLER, FIX, E I GMAX , I E R , NN ) 

DIMENSION A(l),B{l),R(l),Q(l),V(l),E(l),Pll),F(l), 

1 Dll) ,L(i) ,LL(1),MI(1) 



IF 


l 


JB.LT.l .OR. 


J3.GT.IA 


) 


CO 


TO 


6 


IF 


( 


IER.LT.0 .OR. 


IER.GT.7 


) 


GO 


TO 


6 


IF 


l 


NTRY.LT.2 .OR. 


NTRY. GT. 150 


) 


GO 


TO 


6 


GO 


TO 


15 













6 IER = 6 
NN = 0 
RETURN 
15 CONTINUE 

DIVCK = 0.1E-30 

NN =0 

IZ = IIA+D/2 

I A I = IA + 1 

IA2 = IA + 2 

MM = IA*(IA+l)/2 

IAS = I A* I A 

MMP = MM-MM 

EIGMAX = - ABS(EIGMAX) 

IF ( IA.EQ.l ) GO TO 805 
IF l IER.GT.3 ) GO TO 375 
C 

C IF CONTROLLABILITY CHECK REQUESTED, FORM AUGMENTED 

C CONTROLLABILITY MATRIX AND DETERMINE ITS RANK. 

KL = I A* JB 
DO 300 I = 1 , KL 
Dll) = B ( I ) 

300 P(I) = BlI) 

IR = KL 
DO 320 M = 2 , I A 
LI = 0 
LJ = - IA 
DO 310 J= 1 , J B 
LJ = LJ + IA 
DO 310 I =1,1 A 
KI = I - IA 
KJ = LJ 
LI = LI + 1 
FILI) = O.OEO 
DO 310 K=1 , I A 
KI = KI + IA 
KJ = KJ + 1 

310 FILI) = FILI) + A(K1)*P(KJ) 

DO 320 1=1, KL 
IR = IR + 1 
Dl IR) = FI I) 

320 P(I) = F(I) 

CALL GMRANK l D , I A , KL , TOL ER , K ) 

IF l K.EQ.IA ) GO TO 375 
NN = K 
2 IER = 2 
RETURN 

C END CHECK ROUTINE 

C 

375 IF l KQ.EQ.O ) GO TO 35 
DO 20 1=1, MM 
20 Dll) = Q ( I ) 

IR = 0 

KJ = 1 

DO 30 J=l, IA 
KJ = KJ + J - 1 
KI = KJ 

DO 30 I = J , I A 
KI = KI + I - 1 
IR = I R + 1 

30 Q ( I R ) = Dl KI ) 
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C Q = (LOWER TRIANGULAR) Q 

C 

35 IR = 0 

KL = —2*1 A 
KJ = 0 

DO 40 J = 1 , I A 
KL = KL + I A I 
KI = KL 

KJ = KJ + J - 1 
DO 40 1 = J » I A 
IR = IR + 1 
KI = KI + IA 
L { K I ) = IR 

40 L ( I R+K J ) = IR 

IF { JB.GT.l ) GO TO 45 
D{ 1 ) = l.OEO/RU) 

GO TO 55 
45 KL = J 6* JB 
DO 50 1=1, KL 
50 D ( I ) = R ( I ) 

CALL MINV (D, JB,DF,LL,MI) 

55 IR = 0 

DO 60 J = l, IA 
LJ = J - IA 
DO 60 I = J , I A 
LI = I - IA 
KJ = LJ 
IR = IR + 1 
E( IR) = O.OEO 
KL = 0 

DO 60 M = 1 , J 8 
KI = LI 
KJ = KJ + IA 
DO 60 K=1 , JB 
KI = KI + IA 
KL = KL + 1 

60 E { I R ) = E ( I R ) - B(KI)*D(KL)*B(KJ)*1.0D0 
C E = (LOWER TRIANGULAR) - B*R { INVERSE) *B • 

C 

IR = 0 

DO 70 J=l, IA 
LJ = J - IA 
DO 7 0 I = J , I A 
KI = I - IA 
KJ = LJ 
IR = IR + 1 
FUR) = O.OEO 
DO 70 K = 1 , J 3 
KI = KI + IA 
KJ = KJ + IA 

70 FUR) = F ( IR ) + B(KI )*B(KJ)*l.ODO 
DO 80 1=1, MM 
80 F(I) = 2 . OE 0*F { I ) 

IR = 0 

Vll) = O.OEO 
SAVE = O.OEO 
DO 90 J=1 , IA 
KI = J - IA 

IF { VU).LT.SAVE ) V(l) = SAVE 

SAVE = O.OEO 

DO 90 1=1, IA 

KI = KI + IA 

IR = IR + 1 

PUR) = A(KI ) 

90 SAVE = SAVE + ABS ( A ( K I ) ) 

IF { V(1).LT.SAVE ) Vll) = SAVE 
SAVE = FIX*V<1)* S QRT ( FLOAT ( I A ) ) 

C 

C SAVE = FIX^SQRTI I A ) U MAX OVER I OF SUM OVER J ABS(A(I,J))) 
C 
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KI = — I A 
DO 100 1=1, 1A 
KI = KI + I A I 

100 P t K I ) = P ( KI ) + SAVE 
C SOLVE 

C (A+SAVE*I )*S + S*( A+ SAVE*I ) ' = 2*B*B' 

C INITIAL GUESS V { 0 ) = B' *S { INVERSE ) 

C 

CALL MLIAPS ( P , D , L , F , I A , MM , I S , MMP , I A2 ) 

IF ( IS.NE.O ) NN = 1 
DO 110 1 = 1, IAS 
110 P( I) = F(L( I ) ) 

CALL MINV ( P , I A , DF , LL ,M I ) 

IR = IAS 
LJ = - IA 
DO 120 J = 1 , I A 
LJ = LJ + I A 
KI = 0 

DO 120 1=1, JB 
KJ = LJ 
IR = IR + 1 
D( IR) = O.OEO 
DO 120 K = 1 , 1 A 
KI = KI + 1 
KJ = KJ + 1 

120 D( I R ) = D ( I R ) + BCKI)*P(KJ)*1.0D0 
IR = 0 

LJ = IAS - JB 
DO 130 J = 1 , 1 A 
LJ = LJ + JB 
DO 130 1 = 1, IA 
KI = I - IA 
KJ = LJ 
IR = IR + 1 
P ( I R ) = A ( I R ) 

DO 130 K= 1 , J B 
KI = KI + I A 
KJ = KJ + 1 

130 P ( I R ) = P ( IR J - B{ K I )*Q(KJ )*1. ODO 

IF { IER.CT.l .AND. IER.L1.4 .OR. IER.GT.5 ) GO TO 525 
C 

C IF INITIAL GUESS STABILITY CHEfC K REQUESTED, FORM 

C CONTROL MATRIX AND CHECK ITS EIGENVALUES. 

KI = - IA 
V ( 2 ) = O.OEO 
DO 500 1 = 1 , IA 
KI = KI + I A I 
500 V { 2 ) = V ( 2 ) + P ( KT ) 

IF ( V ( 2 ) .GT. IA*EIGMAX ) GO TO 3 
DO 510 1=1, IAS 
510 D ( I ) = P ( I 1 

CALL HSBG (IA,D,IA) 

CALL ATEIG ( I A , D , V , F , LL , I A J 
DO 520 1=1, IA 

IF ( V( I) .GT.EIGMAX ) GO TO 3 
520 CONTINUE 
GO TO 525 
3 IER = 3 
KL = 3 
GO TO 275 

C END CHECK ROUTINE 

C 

525 CONTINUE 

DO 140 1=1, MM 
140 Q( I ) = -Q( I) 
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IR = 0 

LJ = IAS - JB 

DO 150 J = 1 , I A 

LJ = LJ + J6 

DO 150 1 = J , I A 

LI = (I-1)*JB + IAS 

KJ = LJ 

IR = IR + 1 

F( IR) = Q( IR) 

KL = 0 

DO 150 M = 1 , J B 
KI = LI 
KJ = KJ + 1 
DO 150 K=1,JB 
KI = KI + 1 
KL = KL + 1 

150 F ( I R ) = F ( I R ) - D(KI)*R(KL)*D(KJ)*1.0D0 

C SOLVE 

C (A-B*V(0) ) **V( 1) + V(l)*(A-B*VO) ) = - Q - V(0)*R*V(0) 

C 

CALL ML I APS ( P , D, L , F , I A , MM , I S , MMP, I A2 ) 

IF ( IS.NE.O ) NN = NN + 1 
DO 160 1=1, MM 

IF ( F ( I ) .EQ.O.OEO ) F ( I ) = DIVCK 

IF ( ABS(F( I) ) .LT.DIVCK ) F ( I } = S IGN (DIVCK, F( I } } 

160 V( I ) = F ( I ) 

DO 255 M=2,MTRY 
KL = 0 
IR = 0 
LJ = - IA 
DO 210 J = 1 , 1 A 
LJ = LJ + IA 
KI = 0 

DO 210 I = 1,1 A 
KJ = LJ 
IR = IR + 1 
PUR) = O.OEO 
DO 210 K= i , I A 
KI = KI + 1 
KJ = KJ + 1 

210 PUR) = PUR) + E< L(KI ) )*V( L(KJ) )*1.0D0 
IR = 0 
LJ = - IA 
DO 220 J = 1 , I A 
LJ = LJ + IA 
DO 220 I = J , I A 
KI = ( I — 1 ) * I A 
KJ = LJ 
IR = IR + 1 
Ft IR I = Q( IR) 

DO 220 K = 1 , I A 
KJ = KJ + 1 
KI = KI + 1 

220 FUR) = FUR) + V ( L ( K I ) ) *P ( KJ ) *1 .0 DO 
DO 230 1=1, IAS 
230 P ( I ) = A(I ) + P ( I ) 

C SOLVE 

C ( A+E*V(M) ) '*V(M+1) + V(M+1)*(A+E*V(M)} = - Q + V(M)*E*V<M) 

C 

CALL MLIAPS ( P , D , L , F , I A , MM , I S , MM P, I A2 ) 

IF ( IS.EQ.O ) GO TO 235 
NN = NN + 1 

IF ( NN.GT.IZ ) GO TO 7 
235 DO 240 1=1, MM 

IF ( F ( I ) .EQ.O.OEO ) F ( I ) = DIVCK 

IF { ABSIFt I )) .LT.DIVCK ) F(I) = S IGN ( 0 I VCK , F ( I ) ) 

IF { ABSU.OEO - F( I ) /V( I ) ) .GT.TOLER ) GO TO 245 
240 V ( I ) = F ( I ) 

GO TO 265 
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245 DO 250 J = I , MM 

IF ( F(J) . EQ. 0. 0E0 I F ( J ) = DI VCK 

IF ( ABS(FIJ) ) .LT.DIVCK } F(J) = S IGN( D I VCK , F ( J ) ) 
250 V(J) = F ( J } 

255 CONTINUE 
KL = 1 

265 CONTINUE 

DO 270 1=1, MM 
270 Q ( I ) = -Q ( I ) 

275 IF ( KO.EQ.O } GO TO 295 
DO 280 1=1, MM 
280 D( I } = Q ( I ) 

IR = 0 
KJ = 1 

DO 290 J = 1 , 1 A 
KJ = KJ + J - 1 
K I = KJ 
DO 290 I=J , IA 
KI = KI + I - 1 
IR = IR + 1 
290 Q ( K I ) = D( IR) 

295 CONTINUE 

DO 999 1=1, IAS 
999 P(I) = V ( L ( I ) ) 

IF ( I ER-2*( IER/2) .NE.O ) GO TO 645 
C 

C IF POSITIVE-DEFINITE CHECK REQUESTED, ATTEMPT TO 

C FACTOR THE SOLUTION BY CHOLESKY SQUARE ROOT METHOD. 
IF ( P ( 1 } .LT.DI VCK } GO TO 4 
IF ( IA.GT.2 ) GO TO 605 

IF ( P(4)-P(3)*P(3)/P<1) .LT.DIVCK ) GO TO 4 
GO TO 645 
605 IR = 1 

SAVE = SQRTIP(l)) 

DO 610 J = 2 , I A 
IR = IR + IA 

610 D( IR) = P( IR)/SAVE 
IR = 1 
I = 2 

615 IR = IR + I A I 
SAVE = P ( IR) 

LI = IR - I 
DO 620 K=2 , 1 
LI = LI + 1 

620 SAVE = SAVE - D(LI)*D(LI) 

IF ( SAVE. LT.DIVCK ) GO TO 4 
IF ( I .GE . IA ) GO TO 645 
SAVE = SQRT (SAVE) 

LI = IR - I 
LJ = LI 
12 = IR 
M =1 + 1 
DO 640 J=M , I A 
KI = L I 
LJ = LJ + IA 
KJ = LJ 
IZ = IZ + IA 
D ( I Z ) = P ( IZ) 

DO 630 K=2 , I 
KI = KI + 1 
KJ = KJ + 1 

630 D ( I Z ) = D(IZ) - D(KI )*D(KJ)*1.0D0 
640 D ( I Z ) = DIIZ)/SAVE 
1 = 1 + 1 
GO TO 615 
4 IER = 4 + KL 
RETURN 

C END CHECK ROUT INE 

C 
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on 



645 IER = KL 
RETURN 
7 IER = 7 
KL = 7 
GO TO 265 

FOR IA = 1, SOLVE THE RESULTING QUADRATIC EQUATION. 

805 KL = 0 

IF ( IER.LT.4 .AND. ABS ( B { 1 ) ) . LT. DI VCK ) GO TO 2 
D { 1 ) = R(I)/B(1 }*B(1 ) 

D( 2) = A( 1)*D( 1) 

P(l) = D( 2 ) + SQRTI D(2)*D(2) + Q(l)/D(l) ) 

IF ( IER-2*( IER/2) .EQ.O .AND. P(l) .LT.0.0E0 ) GO TO 4 

IER = KL 

RETURN 

END 



SUBROUTINE MLIAPS ( P , D , L , F , I A , MM , I S , MMP, I A2 ) 
DIMENSION P( 1) ,0( 1) ,L( 1) ,F( 1) 

DO 5 1=1, MMP 
5 D(I) = 0.0 EO 

IR = 0 

DO 10 1=1 , IA 
LI = I - IA 
DO 10 J=l, IA 
KJ = J - IA 
KI = LI 
IR = IR + 1 
DO 10 K=1 , IA 
KI = KI + IA 
KJ = KJ + IA 

LJ = L ( K I ) + (L(KJ)-1)*MM 
10 D ( LJ) = D(LJ) + P( IR) 

KI = - IA 
DO 15 1 = 1 , IA 
KI = KI + IA2 - I 
KJ = KI - MM 
DO 15 J = 1 » MM 
KJ = KJ + MM 

15 D ( K J ) = 2 . OE 0*D { K J ) 

CALL SIMQ ( D , F , M M , IS) 

RETURN 

END 
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OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 



SUBROUTINE GMRANK 



PURPOSE 

TO DETERMINE THE NUMERICAL 
A GENERAL MATRIX. 



ROW (COLUMN) RANK OF 



USAGE 

CALL 



GMRANK (D, IA,KL, TOLER, K) 



DESCRIPTION 

D 

IA 

KL 

TOLER - 
K 

REMARKS 

NONE 

SUBROUTINES 

NONE 



OF PARAMETERS 

GENERAL INPUT MATRIX (DESTROYED). 
ROW DIMENSION OF MATRIX D. 

COLUMN DIMENSION OF MATRIX D. 

INPUT CONSTANT USED AS A RELATIVE 
TOLERANCE FOR LOSS OF SIGNIFICANCE. 
ON RETURN, K = RANK OF MATRIX D. 



AND FUNCTION SUBPROGRAMS REQUIRED 



REFERENCES 

PENNINGTON, RALPH H., "INTRODUCTORY COMPUTER 
METHODS AND NUMERICAL ANALYSIS", PUBLISHED BY 
MACMILLAN COMPANY, NEW YORK, 1965. 



SUBROUTINE GMRANK ( D , I A , KL , TOL ER , K ) 
DIMENSION D ( 1 ) 

I A I = I A + I 
NN = MINCH IA,KL) 

IR = - IA 
K = 1 

325 CONTINUE 

IR = IR + I A I 
SAVE = ABS(DUR)) 

LI = K 
LJ = K 
KJ = K - 1 
KI = IR - K 
DO 330 J=K ,KL 
KI = KI + KJ 
DO 330 I =K , I A 
KI = KI + 1 

IF ( ABS(D(KI ) ) .LE. SAVE ) GO TO 330 
LI = I 
LJ = J 

SAVE = ABS ( D { K I ) ) 

330 CONTINUE 

IF ( SAVE.LT .0.1E-35 ) GO TO 365 

IF ( K.GE.NN ) GO TO 375 

IF ( LI.LE.K ) GO TO 345 

KI = IR - IA 

KJ = KI + LI - K 

DO 340 J=K , KL 

KJ - KJ + IA 

KI = KI + IA 

SAVE = D( KI ) 

D ( K I ) = D(KJ) 

340 D ( K J ) = SAVE 

345 CONTINUE 

IF ( LJ.LE.K ) GO TO 355 
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K I = I R — 1 

K J = K I + ( L J-K ) *1 A 

DO 350 I =K » I A 

KI = K I + 1 

KJ = KJ + 1 

SAVE = D ( KI ) 

D{ KI ) = D{ KJ J 
350 D { K J ) = SAVE 

355 CONTINUE 

LI = K + 1 
LJ = IR 

KJ = I R + I A - K 

DO 360 J=LI , KL 

LJ = LJ + IA 

KI = IR 

KJ = KJ + K 

D( LJ ) = D ( L J ) / D ( IR) 

DO 360 I =L 1 » IA 
KI = KI + 1 
KJ = KJ + 1 
SAVE = D( LJ)*D(KI ) 

D( KJ } = D ( KJ ) - SAVE 

IF < ABS(D(KJ ) ) .GE. TOLER* ABS(SAVE) ) GO TO 360 
D ( K J ) = O.OEO 

360 CONTINUE 
K = K + 1 
GO TO 325 
365 CONTINUE 

K _ K 1 

375 CONTINUE 
RETURN 
END 
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